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Abstract 

A  test  bed  is  described  for  the  evaluation  of  timescale  algorithms.  Test  results  for 
two  algorithms  (ATI  and  KAS-l)  using  simulated  data  are  presented.  Both  algorithms 
are  optimum  in  the  sense  that,  with  the  appropriate  selection  of  clock  weights,  they 
are  able  to  produce  an  output  with  minimum  squared  time  prediction  error  in  steady 
state  for  ensembles  of  clocks  whose  noises  are  all  proportional  to  the  noise  of  one  of 
the  members.  The  KAS-l  algorithm  is  also  optimum  during  known  transients  such 
as  startup,  utilizes  robust  outlier  detection,  and  has  the  property  of  time  continuity. 

As  a  result  of  these  improvements,  KAS— 1  is  more  suited  to  automated  operation. 
Neither  algorithm  is  optimum  at  all  sampling  times  for  ensembles  of  clocks  having 
radically  different  noise  processes.  The  KAS— 1  algorithm  and  a  clock  simulator,  both 
representing  a  new  Kalman-filter  approach,  are  documented  herein. 

INTRODUCTION 

The  Master  Clock  Upgrade  Program  at  USNO  is  a  SPA  WAR-  and  IJSNO-funded  project  which  tasks 
NRL  and  USNO  to  procure,  monitor,  and  evaluate  advanced  frequency  standards  (trapped  mercury 
ion  devices  and  hydrogen  masers)  currently  under  development  by  others;  devise  and  test  a  low- 
noise  time-difference  measurement  system  for  all  the  cesium  clocks  and  hydrogen  masers  at  USNO; 
and  develop  a  test  bed  for  algorithms  designed  to  generate  an  optimum  timescale^.  In  order  to 
accomplish  the  latter,  a  software  package  is  being  developed  by  NRL  that  will  be  an  integral  part 
of  the  data  collection  and  management  system  at  USNO.  The  current,  version  of  the  software  that 
comprises  the  test  bed  is  written  in  “Rocky  Mountain  Basic”  for  the  IIP  (Hewlett  Packard)  series  300 
computer.  The  next  version  of  the  test  bed  package  will  be  in  the  standard  “C”  language  and  will  be 
less  machine  dependent.  In  addition  to  the  standard  clock  analysis  tools,  the  test  bed  will  allow  two 
algorithms  to  run  simultaneously  using  the  same  data,  which  may  be  real  or  simulated.  Tn  the  case  of 
the  latter,  the  output  of  an  algorithm  is  compared  to  truth,  where  truth  is  known  by  definition. 

As  a  joint  project,  USNO  and  NRL  are  making  complementary  contributions  to  the  algorithm  test  bed. 
NRL’s  contribution  is  the  experience  in  developing  automated  clock  systems.  USNO’s  contribution  is 
the  experience  of  and  knowledge  about  maintaining  a  timescale. 

‘Work  supported  in  part  by  Naval  Research  Laboratory,  Contract  Nos.  NOOl 73-87-M — J 780,  N00173-88-M  0031, 
N00173-88-M-5999,  and  N00173-88-M-G568 
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NRL  is  developing  a  number  of  automated  clock  systems  for  field  use  designed  to  have  the  following 
features: 


1.  the  systems  are  fully  automated  clock  ensembles  with  real-time  output; 

2.  they  can  be  steered  to  UTC  (USNO)  automatically  via  GPS; 

3.  their  timescale  algorithms  permit  the  easy  addition  or  deletion  of  clocks;  and 

4.  the  algorithms  generate  a  timescale  more  uniform  than  any  single  clock  in  the  ensemble. 

For  our  purposes,  we  will  define  a  timescale  as  a  multifaceted  procedure  defining  the  scale  of  time  by 
means  of  an  algorithm  that  usually  entails  filtering  or  smoothing  the  data  from  two  or  more  clocks  or 
clock  ensembles,  and  possibly  steering  one  or  more  of  the  clocks  toward  the  ensemble  or  steering  them 
or  the  ensemble  itself  toward  a  more  stable  reference.  We  define  a  clock  ensemble  to  be  a  procedure  by 
which  a  single  output  is  obtained  from  two  or  more  clocks.  Generally  an  ensemble  performs  better  than 
its  individual  members.  For  a  set  of  N  clocks  with  equal  noise  processes,  an  output  can  be  obtained 
that  is  y/N  more  stable  than  a  single  clock,  assuming  the  absence  of  intercorrelations  and  systematic 
errors.  For  a  set  of  clocks  whose  noises  are  not  proportional,  the  output  is  more  complicated,  as  we 
shall  see. 


RESOURCES 

As  mentioned  above,  the  test  bed  uses  both  simulated  and  real  data.  One  source  of  the  latter  is  the  2~ 
picosecond-noise  time  measurement  system  developed  by  NRL  for  USNO  that  is  connected  to  at  least 
twenty  high-performance  cesium  clocks,  twelve  hydrogen  masers,  and  three  mercury  stored-ion  devices 
(SIDs).  The  measurement  system  incorporates  a  phase-difference  measurement  subsystem  developed 
by  the  first  author  while  at  NBS^l  The  cesiums  and  masers  provide  phase- difference  measurements, 
while  the  SIDs  give  frequency  measurements.  This  mix  of  measurement  types,  a  challenge  to  most 
ensembling  procedures,  is  handled  quite  naturally  by  a  Kalman-filter  approach.  In  addition  to  these 
new  data,  there  is  available  the  old  USNO  database  of  several  years  from  two  masers  and  about  fifty 
cesiums, 

NRL  has  available  for  this  effort  a  database  of  phase  measurements  from  three  masers,  at  least  ten 
(some  high-performance  and  some  GPS-type)  cesiums,  and  several  rubidium  clocks.  Tn  addition,  NRL 
has  a  means  of  measuring  the  USNO  Master  Clock  to  the  50-picosecond  level  via  the  carrier  frequency 
of  a  local  television  station  using  a  system  similar  to  one  devised  in  France^'  1.  This  precise  method 
of  remote  measurement  will  permit  clocks  several  miles  apart  to  be  compared  and  allow  the  potential 
ensembling  of  clocks  at  both  NRL  and  USNO. 

The  clock  simulator  employed  in  this  test  bed  is  a  Kalman-filter  type  introduced  in  the  next  section. 
Simulated  data  have  several  advantages  over  real  data: 

1.  any  amount  of  data  over  any  length  of  time  can  be  generated; 

2.  data  free  of  systematic  errors  or  containing  specific  types  of  errors  can  be  produced; 

3.  any  number  of  clocks  or  combination  of  clock  types  can  be  modelled;  and 

4.  most  importantly,  truth  is  known  by  definition. 

To  illustrate  the  latter,  consider  the  simulation  of  a  single  clock,  which  involves  the  generation  of  a 
series  of  numbers  that  represent  the  accumulation  of  phase,  usually  over  equal  time  intervals.  The 
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numbers  are  produced  with  a  clock  model  in  mind  and  against  which  they  may  be  verified  statistically. 
To  generate  phase-difference  data,  two  or  more  clocks  are  simulated  simultaneously;  one  clock  is  then 
selected  as  the  reference  and  its  phase  value  at  each  measurement  step  is  subtracted  from  the  other 
clocks. 

THE  ALGORITHMS 

The  first  algorithm  tested  was  ATI,  developed  at  NIST  * The  source  code  was  provided  by  NBS 
(now  NIST)  in  1987.  It  is  an  ad  hoc  procedure  designed  for  a  steady  state  ensemble  of  clocks  manifest¬ 
ing  the  types  of  noise  typical  of  cesium  clocks.  Exponential  filters  are  used  in  the  frequency  estimation 
and  weight  calculation. 

The  second  algorithm,  called  the  Kalman  Aiding  Sources  Algorithm  (KAS-1),  developed  in  1987  by 
Ball  Aerospace  for  NRL,  uses  Kalman-filter  methodology.  The  KAS-1  algorithm  was  first  employed  in 
an  automated  clock  ensemble  at  the  Master  Control  Station  of  the  Global  Positioning  System  (GPS) 
at  the  Consolidated  Space  Operations  Center  (CSOC),  Falcon  AFB  and  was  presented  at  the  Third 
International  Time  Scale  Algorithm  Symposium  held  in  Torino  in  1988^1  This  paper  provides  a  more 
complete  description  as  well  as  operational  and  simulated  performance  data.  The  KAS  -1  approach  is 
now  being  tested  by  NIST  for  use  with  ATI  and  their  results  are  described  in  the  paper  by  Weiss  and 
Weissert  in  these  same  proceedings.  The  KAS-1  algorithm  was  the  first  successful  use  of  a  Kalman 
filter  in  connection  with  clock  ensembles,  and  this  approach  has  been  enhanced  for  use  at  USNO. 

KALMAN  APPROACH  TO  SIMULATING  AND  ENSEMBLING 
CLOCKS 

Background 

The  ensembling  method  employed  in  KAS— 1  uses  one  clock  to  form  ari  initial  estimate  of  time  and  uses 
each  additional  clock  to  improve  or  ‘aid’  the  initial  estimate.  The  same  approach  is  used  in  the  NIST 
ATI  timescale.  The  work  described  here  improves  and  extends  the  ATI  algorithm  in  several  areas: 
optimum  filtering  of  all  clock  types,  automatic  startup  and  clock  addition,  robust  outlier  processing, 
time  and  frequency  step  detection,  real  time  clock  parameter  estimation,  and  adaptive  modelling. 

The  KAS-1  algorithm  utilizes  Kalman  filters  to  provide  the  necessary  state  estimation  and  forecasting 
functions,  resulting  in  estimates  which  are  optimum  in  the  minimum  squared-error  sense  both  in  steady 
state  and  during  transient  conditions  such  as  turn-on.  Any  mix  of  clocks  can  be  included  provided  that 
the  correct  models  and  covariances  are  used.  The  use  of  the  Kalman  filter  is  a  computational  device 
and  any  other  method  of  minimum  squared-error  state  estimation  based  on  the  same  assumptions 
and  system  dynamics  would  produce  an  identical  result. 

In  contrast,  the  NIST  ATI  algorithm  uses  fixed  length  exponential  filters  to  perform  frequency  es¬ 
timation.  It  is  optimum  only  in  steady  state  and  requires  external  clock  calibration  for  startup  and 
addition  of  clocks  to  the  ensemble. 

Definition  of  the  KAS— 1  Ensemble 

Each  of  the  N  clocks  provides  an  independent  forecast  of  the  time  of  clock  r,  the  reference  clock,  with 
respect  to  the  ensemble  at  time  t  +  8.  These  individual  forecasts  are  called  ‘primitive’  forecasts  of  the 
reference  clock  time.  The  primitive  forecast  of  the  time  of  the  reference  clock  at  t  +  S  using  clock  i  is 

ure{t  T  ^  “b  “  uie(t  +  8\t)  -f-  Ur{  (t  \  8 )  (l) 
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where  t qe(t  8 |t)  is  the  forecast  of  the  time  of  clock  i  with  respect  to  the  ensemble  at  time  t  f  8  based 
upon  the  true  state  of  clock  t  with  respect  to  the  ensemble  at  time  t  and  ur,(t  +  8)  is  the  clock  r  vs. 
clock  i  time  difference. 

The  KAS-1  algorithm  utilizes  the  natural  definition 

N 

Ure(t  +  8)  =  y]  a,(t  +  8)u're(t  +  8\t  +  8).  (2) 

»'=i 

The  weights,  a,(t),  may  be  chosen  in  any  way  subject  only  to  the  restriction  that  the  sum  of  the  weights 
is  one.  In  the  absence  of  noise,  the  ensemble  definition  is  deterministic.  However,  in  the  presence 
of  clock  and  measurement  noise,  the  computational  problem  is  to  find  the  minimum  squared-error 
estimate  of  the  time  of  the  reference  clock  with  respect  to  the  ensemble. 

Computational  Methodology  for  the  KAS-1  Algorithm 

One  method  of  estimating  ensemble  time  is  to  use  the  definition  of  Eq.  2  directly.  Taking  the  estimate 
of  both  sides,  one  obtains 

N 

ure(t  +  £|f  +  8)  =  a,(t  +  6)[u,e(t  +  8\t)  +  ur,-(t  +  <5)].  (3) 

t=i 

The  estimate  of  the  time  of  clock  i  with  respect  to  the  ensemble  at  time  t  +  8  based  on  measurements 
through  time  t,  the  first  term  in  the  square  brackets  on  the  right  side  of  Eq.  3,  is  calculated  from  the 
results  of  the  last  computation  of  the  ensemble.  Thus 

8 2 

uie(t  +  <5|t)  =  z,e(t|t)  +  8yie(t\t)  +  ~wie(t\t)  (4) 

where  x  is  the  phase  before  perturbation  by  white  phase  noise,  y  is  the  frequency  and  w  is  the  frequency 
aging.  The  second  term  in  the  square  brackets  on  the  right  side  of  Eq.  3  is  the  minimum  square-error 
estimate  of  the  clock  pair  difference  state.  This  estimate  is  easily  computed  using  the  appropriate 
Kalman  filter. 

Kalman  Filter  for  Computing  Clock  Difference  States 

The  n  parameters  which  are  to  be  estimated  are  formed  into  an  n  dimensional  state  vector  x(t).  The 
system  evolves  from  time  t  to  time  t  +  8  according  to 

x(t  +  8)  =  *(«)*(*)  +  TS(t  +  8\t)  +  $(<5)p(i),  (5) 

where  the  nXn  dimensional  state  transition  matrix  <&(£)  embodies  the  system  model,  the  n  dimensional 
vector  Ts(t  +  £|t)  contains  the  noise  inputs  to  the  system  during  the  interval  from  t  to  t  +  5,  and  the 
n  dimensional  vector  p(t)  contains  the  control  inputs  made  at  time  t.  The  state  transition  matrix  is 
assumed  to  depend  on  the  length  of  the  interval,  but  not  on  the  origin.  Each  element  of  s(t  +  8\t) 
is  Normally  distributed  with  zero  mean  and  is  uncorrelated  in  time,  thereby  generating  a  random 
walk  in  the  elements  of  the  state  vector.  The  observation  vector  z(l)  is  described  by  the  measurement 
equation 

z(t)  =  H(f)x(t)  +  v(t).  (6) 

The  r  observations  made  at  time  t  are  related  linearly  to  the  n  elements  of  the  state  vector  by  the  r  x  n 
dimensional  measurement  matrix  H(t)  and  the  r  dimensional  white  noise  vector  v(i).  Kalman  and 
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Bucy  defined  a  recursive  procedure  for  estimating  the  next  state,  which  requires  the  mean  squared- 
error  of  the  estimates  from  the  true  state  to  be  minimum. 

The  error  in  the  estimate  of  the  state  vector  after  the  measurement  at  time  t  is  x(z|i) - x(i)  and  the 

error  covariance  matrix  is  defined  to  be 

P(t|f)  =  E  |  [x(t|t)  -  -x(t)]  [x(tjz)  -  -x(f)]  J  .  (7) 

The  diagonal  elements  of  this  n  x  n  matrix  are  the  variances  of  the  estimates  of  the  components  of 
x(f)  after  the  measurement  at  time  t.  Next,  the  error  covariance  matrix  just  prior  to  the  measurement 
at  time  t  +  8  is  defined  as 

P(Z  +  £|Z)  =  £'|jx(Z  +  ^|Z)  -x(z  +  £)j  [x'(Z  +  <5jz)  -  x(t- j-  <5)j  |  .  (8) 

Finally,  R(Z)  is  the  covariance  matrix  of  the  measurement  noise  and  Q(Z  +  5|Z)  is  the  covariance  matrix 
of  the  system  noise  generated  during  the  interval  from  t  to  t  +  8 

R(Z)  =  E  [v(Z)v(/.)T]  (9) 

Q (t  +  8\t)  =  E  [s(t  +  §|i)s(Z  +  S\t)T]  .  (10) 

The  error  covariance  matrix  evolves  according  to  the  system  model 

P (t  +  8\t)  -  &(8)?(t\t)$(8)T  +  TQ(Z  +  <5j*)rT.  (11) 

The  new  estimate  of  the  state  vector  depends  on  the  previous  estimate  and  the  current  measurement 

x(Z  +  £|f  +  6)  =  #(<5)x(Z|Z)  +  $  (6)p(Z) 

+  K {t  +  8)  [z(Z  +  5) -H(Z  +  5)$(^(Z|Z) -H(Z  +  5)$(^)p(f)]  ,  (12) 

where  the  gain  matrix,  K(t  +  8 ),  determines  how  heavily  the  new  measurements  are  weighted.  The 
optimum  or  Kalman  gain,  K opt,  is  determined  by  minimizing  the  “square  of  the  length  of  the  error 
vector,”  i.e.,  the  sum  of  the  diagonal  elements  (the  trace)  of  the  error  covariance  matrix 

Kopt(t  +  8)  =P{t-\-S\t)U{t  +  8)T  [H(t+£)P(Z  +  <5|f)H(i  i-  8)T  +  R(i  +  5)] _1 .  (13) 

Finally,  the  updated  error  covariance  matrix  is 

P(t  +  8lt  +  8)  =  [I-K(t  +  8)H(t  +  8)]P(t  I  8jt)[I-K(/,  +  8)H(t  t  8)]r 

+  K(t  +  8)R(t  l-8)K(t  +  8)T  (14) 

where  I  is  the  identity  matrix.  If  the  optimum  filter  gain  is  used,  Eq.  1 4  reduces  to  a  simpler  form: 

P (t  -I-  8\t  +  8)  •  [I  -  K opt(t  t  8)U(t  +  8)}  P(i  j  S\t).  (15) 

Equations  11  through  15  define  the  Kalman  filter,  and  so  defined  it  is  an  optimal  estimator  in  the 
minimum  squared  error  sense.  Each  application  of  the  Kalman  recursion  yields  an  estimate  of  the  state 
of  the  system  which  is  a  function  of  the  elapsed  time  since  the  last  filter  update.  Updates  may  occur  at 
any  time.  In  the  absence  of  observations,  the  updates  are  called  forecasts.  Simultaneous  observations 
may  be  processed  in  parallel  or  serially,  in  any  order,  if  the  measurement  noise  is  uncorrelated  with 
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the  process  noise.  The  interval  between  updates,  8,  is  arbitrary  and  is  specifically  not  assumed  to  be 
constant. 

The  estimates  of  the  clock  states  relative  to  clock  r  are  obtained  from  N  —  1  independent  Kalman 
filters.  The  four  dimensional  state  vectors  are: 

Uri  (t) 

mI!  ■  <“> 

wri  (t)  _ 

Every  clock  pair  has  the  same  state  transition  matrix 


015  82/2  ' 
0  1  8  82/ 2 
0  0  1  8 
0  0  0  1 


(17) 


and  F  matrix 

1  1  0  0  ' 
0  10  0 
0  0  10 
0  0  0  1 

The  system  covariance  matrices  are: 

qir(t  +  8\t)  = 


(18) 

(19) 


"  s"'(t)fh  0  0  0 

0  5‘r(t)5  +  5;r(i)53/3  +  5;r(t)55/20  Sir(t)82/2  +  S"(t)8Jl/8  S’r(t)53/6 

0  Slr(t)82/2  +  Sir{t)8A/8  S';{t)5  +  Sir{t)8y  3  S"(t)S2/ 2 

0  S';r{t)8s/ 6  Sjr(t)82/ 2  S'(r\t)8 

where  the  clock  pair  spectral  densities  are  the  sum  of  the  individual  contributions  from  each  of  the 
clocks: 

sir  =  S'  +  Sr.  (20) 

The  white  phase  measurement  noise  is  given  by  the  measurement  model: 

zrt  =  Hxn  +  vn ,  (21) 

where  each  time  measurement  is  described  by  the  same  4  x  1  row  matrix 

Hr,  =  (  1  0  0  0  ).  (22) 

The  updated  difference  states  are  given  by  Eq.  12,  but  only  the  time  state  is  used  in  the  following 
steps  of  the  ensemble  calculation. 

Calculation  of  the  Weighted  Average  of  the  Times  of  the  Clocks 

Although  it  is  not  necessary,  a  Kalman  filter  is  used  to  calculate  the  ‘weighted  average’  of  the  N 
primitive  estimates  of  the  time  of  clock  r.  Kalman  filters  arc  normally  applied  to  dynamic  systems, 
but  they  are  equally  relevant  to  the  static  problem  of  ‘averaging’  the  N  estimates  of  the  time  of  clock 
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r  at  sample  time  t.  The  choice  of  the  Kalman  filter  methodology  is  motivated  by  the  fact  that  it 
deals  with  the  clock  weights  in  an  easy  and  natural  manner  and  automatically  provides  scale  factors 
for  use  in  outlier  deweighting.  The  deweighting  scheme  is  similar  to  one  developed  by  Percival^  for 
maximum  likelihood  estimation,  but  it  has  been  adapted  for  use  with  the  Kalman  methodology.  It  is 
desirable  for  the  resulting  estimator  (with  outlier  deweighting)  to  be  both  efficient  and  robust.  For 
certain  types  of  outliers,  the  method  used  here  is  both  more  efficient  and  more  robust  than  rejection 
rule  methods. 

The  state  transition  matrix  is  a  scalar  and  is  equal  to  one  since  the  system  does  not  evolve  between 
‘observations’.  The  ‘measurement’  model  is 

Ure{t)  =  «re(t|t)  +  Vt(t)  (23) 

where  the  random  shocks  v,-(t)  have  variances  a ?(f)  and  arc  used  to  set  the  weights  of  the  clocks.  For 
example,  v,-(t)  =  v(t)  independent  of  the  clock  number  produces  equal  weighting.  On  the  other  hand, 
setting  Uj(t)  equal  to  the  estimated  time  prediction  error  of  clock  i  since  the  last  ensemble  calculation 
produces  an  ensemble  with  minimum  time  prediction  error.  Note  that  the  latter  operating  condition 
is  not  necessarily  the  best  since  the  clocks  may  not  be  perfectly  modeled. 

It  is  possible  to  solve  for  the  Kalman  gain  in  closed  form: 

(24) 

However,  the  full  Kalman  recursion  is  more  suited  to  the  deweighting  of  outliers.  Since  a  scalar  is 
being  estimated,  «&  =  H  =  1,  Q  =  0,  and  E  —  cr|(t).  In  order  to  deweight  outliers,  a  robust  initial 
guess  is  required  for  the  time  of  clock  r  with  respect  to  the  ensemble  after  the  measurement  at  time 
t  +  S.  The  median1  of  the  available  estimates  is  a  good  choice  since  it  is  not  unduly  affected  by  outliers. 

uTe(t  +  S\t  +  £)i  =  median  j u\e(t  -f  S\t  +  £)}  (25) 

This  estimate  is  subsequently  refined  by  using  the  information  in  the  remaining  primitive  estimates. 
The  primitive  estimates  are  sorted  according  to  their  deviation  from  the  median  so  that  clock  I(k ) 
is  the  kth  nearest  the  median  and  /(l)  is  the  median.  To  start  the  Kalman  recursion  one  needs  to 
initialize  the  state  and  the  covariance  matrices.  Having  just  completed  the  estimate  of  the  ensemble 
time  for  the  measurements  taken  at  time  t,  the  ensemble  calculation  for  time  t  +  S  commences  by 
setting  the  initial  estimate  of  the  time  of  clock  r  vs.  the  ensemble  equal  to  the  ‘primitive  estimate’ 
using  clock  /( 1).  The  starting  point  is  important  only  for  outlier  deweighting.  A  ‘robust’  initial 
estimate  for  the  time  of  clock  r  is  required  and  using  one  of  the  primitive  estimates  avoids  biasing  the 
final  estimate.  The  first  estimate  must  be  robust  because  it,  cannot  be  deweighted.  This  implies  no 
loss  of  generality  since  it  is  not  possible  to  detect  an  outlier  with  only  one  clock  and  it,  is  not  possible 
to  assign  responsibility  for  an  outlier  time  difference  measurement  with  only  two  clocks.  That  is,  there 
is  one  more  cycle  in  the  Kalman  recursion  than  is  needed  for  outlier  deweighting.  The  recursion  begins 
by  setting  the  initial  state  estimate: 

ure[t  +  S\t  +  S')  i  =  (/,  +  S\t  +  S') . 


m  = 

1  -jlF 


1When  the  number  of  points  is  odd,  the  number  of  observations  larger  than  the  median  equals  the  number  of 
observations  smaller  than  the  median.  When  the  number  of  observations  is  even,  the  median  is  one  of  the  two  central 
observations.  Other  definitions  of  the  median  are  sometimes  used. 
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The  covariance  matrix  prior  to  processing  the  measurement  of  clock  1(2)  is 

P(/(2)|J(1))  =  ^(1) 

and  the  optimum  Kalman  gain  is 


K 


pm\m 

m  P(/(2)|J(1))  +  *| 


J(2) 


At  this  point  in  the  calculation,  the  standard  Kalman  methodology  must  be  modified  to  account 
for  the  possible  presence  of  outlier  observations.  If  primitive  estimate  1(2)  does  not  come  from  the 
distribution  included  in  the  physical  model,  then  it  should  not  be  used  to  determine  the  Kalman  state 
estimate.  In  order  to  preserve  the  continuity  of  the  state  estimates,  measurements  should  not  be 
rejected  precipitously,  rather  the  Kalman  gain  should  be  reduced  as  the  measurement  departs  from 
expectation.  A  non-optimum  Kalman  gain  is  calculated  from 


,  _  Vtm[?/( 2)] 

A/(2)  -  “  K 


9/(2) 


1(2) 


where 


^HA(q)  =  < 


q 

0 


kl  < a 


if 

if  a  <  \q\  <  b 
if  k|  >  b 


(26) 


(27) 


is  Hampel’s  ip  function  and 


?/(2) 


urg  ^ (t  +  <$|£  +  (5)  —  ure(t  +  S\t  +  ^)i 


S  2 


(28) 


The  scale  factor,  s2,  is  discussed  later.  The  second  step  in  the  Kalman  recursion  is  completed  by  using 
the  modified  Kalman  gain  to  calculate  the  updated  covariance  matrix  according  to  the  formula  for 
non-optimum  Kalman  gain 


P[/(2)|J(2)]  =  [l  —  A}(2)]2  P[J(2)|/(1)]  +  [a;(2)]2^(2) 
and  to  calculate  the  new  estimate  of  the  time  of  clock  1  with  respect  to  the  ensemble 
ure(t  +  <5|£  +  8)2  —  ^1  —  «re(t  +  <5]£  +  5)i  +  R'j^u1^ (t  +  8\t  +  8). 


(29) 


(30) 


The  recursion  continues  in  this  way  until  all  the  measurements  are  processed.  The  computation  for 
clock  I(j)  starts  by  setting 

P[/(i)|/(i-l)]  =  P[/(i-l)|/(j-l)]. 


The  first  time  the  recursion  is  performed,  each  of  the  constants  Sj  is  set  equal  to  the  maximum  of  the  <7y. 
This  prevents  the  high  noise  level  of  the  median  from  causing  a  very  quiet  clock  to  be  downweighted. 
Next,  the  computation  is  repeated  setting  the  initial  estimate  of  the  time  of  the  reference  with  respect 
to  the  ensemble  equal  to  the  primitive  estimate  closest  to  the  result  of  the  previous  computation  and 
setting  Sj  equal  to  the  second  largest  ay  or  its  own  prediction  error,  whichever  is  larger.  The  recursion 
is  repeated  until  each  estimate  is  processed  using  the  <7,-  of  the  corresponding  clock.  If  any  of  the 
clocks  were  downweighted  in  this  process,  then  the  new  estimate  of  the  time  of  clock  r  with  respect 
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to  the  ensemble  is  compared  with  the  previous  estimate  and  the  process  is  repeated  until  the  change 
in  the  estimate  between  iterations  is  negligible.  When  there  are  no  outliers,  this  recursive  solution  is 
identical  to  Eq.  3.  In  the  presence  of  outliers,  it  is  still  a  weighted  average  oT  the  estimates  from  the 
individual  clocks.  The  weights  of  the  clocks  are  given  by 

aHf)  =  K'i(j)  II  U  “  K'i(i)) 

*=j+i 

where  Tfj ^  =  1.  To  preserve  the  reliability  of  the  ensemble,  one  usually  limits  the  weights  of  each  of 
the  clocks  to  some  maximum  value,  amax.  A  deweighting  factor,  k,  may  be  calculated  to  accomplish 
this  goal  by  starting  at  index  I(N )  and  proceeding  to  index  7(1)  so  that  kK1  are  new  weights  which 
do  not  exceed  the  limit.  Each  time  the  Kalman  recursion  is  used  to  estimate  the  time  of  clock  r,  the 
previously  calculated  /c  is  used.  The  process  rapidly  converges  to  a  stable  value. 

Finally,  the  times  of  the  remaining  clocks  are  computed  from  the  time  of  clock  r  and  the  estimated 
time  differences.  Thus 


(31) 


flic  (f  +  6  [  t  +  6 )  —  ure  (t  +  6\t  +  £)  —  uri  (f-j-<5i|£  +  £). 


(32) 


Calculation  of  the  Frequency  and  Frequency  Aging  States 

The  time  of  clock  i  relative  to  the  ensemble  is  used  as  input  to  a  Kalman  filter  that  estimates  the 
frequency  and  frequency  aging  of  clock  t  relative  to  the  ensemble.  This  filter  uses  the  state  transition 
matrix  of  Eq.  17.  It  is  assumed  that  the  ensemble  is  so  large  that  its  contribution  to  the  noise  can  be 
accounted  for  to  first  order  only.  Based  on  this  assumption  of  a  large  ensemble,  the  system  covariance 
matrix  is: 

Qie(t  +  8\t)  =  Qi{t  +  S\l).  (33) 

and 

1  -  <z,(f)  1  -  a,-(t )  0  0 

0  1  -  a,  (t)  0  0 

0  0  10 

0  0  0  1 

The  change  in  F  from  the  case  of  a  single  clock  or  clock  pair  reflects  the  fact  that  each  clock  is  a 
member  of  the  ensemble  with  respect  to  which  it  is  measured. 

The  measurement  model  is  noiseless  since  the  calculated  times  of  the  N  clocks  are  used  as  pseu¬ 
domeasurements.  This  procedure  guarantees  that  the  Kalman  filter  reproduces  the  phase  estimates 
obtained  directly  from  the  ensemble  definition  while  simultaneously  estimating  consistent  values  of 
the  frequency  and  frequency  aging. 

Adaptive  Modelling 

Parameter  Estimation 

A  variance  analysis  technique  compatible  with  irregular  observations  has  been  developed^.  The  vari¬ 
ance  of  the  innovation  sequence  of  the  Kalman  filter  is  analyzed  to  provide  estimates  of  the  parameters 
of  the  filter.  The  result  is 


E 


V (t  +  S)v(t  4-  sf  =  H(£  +  8)F(t  +  £|i)II(i  -l  S)T  +  R(£  +  8). 


(35) 
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Like  the  Allan  variance  analysis,  which  is  performed  on  the  unprocessed  measurements,  the  innovation 
variance  analysis  requires  only  a  limited  memory  of  past  data.  However,  the  forecasts  provided  by  the 
Kalman  filter  allow  the  computation  to  be  performed  at  arbitrary  intervals  once  the  algebraic  form  of 
the  innovation  variance  has  been  calculated.  Adaptive  modelling  begins  with  an  approximate  Kalman 
gain,  K.  As  the  state  estimates  are  computed,  the  variance  of  the  innovations  on  the  left  side  of 
Eq.  35  is  also  computed.  The  right-hand  side  of  this  equation  is  written  in  terms  of  the  actual  filter 
element  values  (covariance  matrix  elements)  and  the  theoretical  parameters.  Finally,  the  equations 
are  inverted  to  produce  improved  estimates  for  the  parameters. 

Using  the  autocovariance  function  to  solve  for  the  parameters  is  inappropriate  here,  because  the  auto¬ 
covariance  function  is  highly  correlated  from  one  lag  to  the  next  and  the  efficiency  of  data  utilization 
is  therefore  small.  Instead,  only  the  autocovariance  of  the  innovations  for  zero  lags,  i.e.,  the  variance 
of  the  innovations,  is  used.  The  variance  is  given  by 

E[uic{t  +  S)i7ie(t  +  6)T]  =  Pi1(t|t)  +  2«P*12(t|t)  +  tf*Pis(*|0 

+  53P*3(t|t)  +  ^P‘22(t|t)  +  ^r*33(t|t) 

+  [Si(t)S  +  4(t)M(l  -  a,(t)]2  +  5*  (i)y  +  4*4  +  <.•(*). 

(36) 

It  is  assumed  that  the  oscillator  model  contains  no  hidden  noise  processes.  This  means  that  each  noise 
in  the  model  is  dominant  over  some  region  of  Fourier  frequency  space.  The  principle  of  parsimony 
encourages  this  approach  to  modelling.  Inspection  of  Eq.  36  leads  to  the  conclusion  that  each  of  the 
parameters  dominates  the  variance  of  the  innovations  in  a  unique  region  of  prediction  time  interval, 
S,  making  it  possible  to  obtain  high-quality  estimates  for  each  of  the  parameters  through  a  process  of 
bootstrapping. 

For  each  parameter  to  be  estimated,  a  Kalman  filter  is  computed  using  a  subset  of  the  data  chosen 
to  maximize  the  number  of  predictions  in  the  interval  for  which  that  parameter  makes  the  dominant 
contribution  to  the  innovations.  The  filters  are  designated  0  through  4,  starting  with  0  for  the  main 
state  estimation  filter  which  runs  as  often  as  possible.  Each  innovation  is  used  to  compute  a  single¬ 
point  estimate  of  the  variance  of  the  innovations  for  the  corresponding  5.  Substituting  the  estimated 
values  of  the  remaining  parameters,  Eq.  36  is  solved  for  the  dominant  parameter,  and  the  estimate  of 
that  parameter  is  updated  in  an  exponential  filter  of  appropriate  length.  If  the  minimum  sampling 
interval  is  too  long,  it  may  not  be  possible  to  estimate  one  or  more  of  the  parameters.  However,  there 
is  no  deleterious  consequence  of  this  situation,  since  a  parameter  that  cannot  be  estimated  is  not 
contributing  appreciably  to  the  prediction  errors.  It  is  not  necessary  to  use  separation  of  variances 
to  estimate  individual  clock  parameters  since  the  T  matrix  takes  into  account  the  first  order  clock- 
ensemble  correlations. 

The  KAS-1  algorithm  is  adaptive  in  two  ways.  First,  the  clock  noise  parameters  are  updated  after  each 
cycle  of  the  recursion  and  used  during  the  next  cycle  of  the  recursion.  Second,  the  clock  weights  are 
updated  after  each  cycle  based  on  the  real  time  parameter  estimates  or  the  measured  time  prediction 
errors.  For  minimum  time  prediction  errors,  the  variances  (inverse  Weights)  needed  for  the  calculation 
of  the  time  of  the  reference  clock  with  respect  to  the  ensemble  are 

<7, 2  =  4  (t)fh  +  Sj(t)«  +  Sl(t)83/ 3  +  Sj(t)Ss/ 20.  (37) 
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Kalman  Clock  Simulation  Algorithm 

The  algorithm  presented  here  may  be  used  to  compute  the  simulated  state  of  a  precision  clock  at  any 
time  in  the  future  based  on  the  current  state  of  the  clock  and  appropriate  noise  inputs.  White  phase 
measurement  noise,  white  phase  additive  noise,  white  frequency  noise,  random  walk  frequency  noise, 
and  random  walk  frequency  aging  noise  may  all  be  included.  Although  the  algorithm  is  quite  simple,  it 
has  several  advantages  over  the  ARIMA  technique  described  by  Barnes^.  The  time  interval  between 
observations  is  included  explicitly  so  that  a  time  series  of  observations  separated  by  varying  intervals 
may  be  created.  In  addition,  the  parameters  are  easily  calculated  from  the  power  spectral  densities  of 
the  noise.  These  physical  parameters  are  the  same  ones  required  for  the  Kalman  analysis  of  the  clock. 
The  final  advantage  is  that  the  simulated  clock  noise  for  each  state  is  appropriately  correlated  with 
the  noise  on  other  states,  as  should  be  the  case  since  the  physical  model  for  a  clock  is  an  integrator. 

The  simulator  is  based  on  the  same  equation  of  motion  for  a  clock,  Eq.  5,  used  in  Kalman  filter 
analysis.  Each  element  of  s {t  +  5|t)  is  Normally  distributed  with  zero  mean  and  is  uncorrelated  in 
time.  The  problem  is  to  generate  simulated  noise  vectors,  s'(f  |-  5\t),  having  the  desired  covariance, 
given  in  Eq.  19.  The  spectral  densities  which  appear  in  the  system  covariance  matrix  may  be  written 
in  terms  of  the  standard  oscillator  noise  coefficients.  The  values  and  units  are  shown  in  Table  1. 

Table  1:  Relationship  between  spectral  densities  and  h  coefficients 


Spectral  density 

h  coefficient 

units 

S0'(t) 

(2x  )~2h2(l) 

seconds3  j 

Se(t) 

MO 

seconds 

Sfj.  (t) 

(2ff)2h-2(t) 

seconds-1 

MO 

(27T ) 4 /l _ 4  ( t ) 

seconds  A 

The  goal  is  to  derive  a  solution  for  the  noise  vector,  s,  in  the  form 

s  =  Ar , 


(38) 


where  r  is  a  four  dimensional  vector  of  Normal  deviates  with  unity  variance  whose  elements  are 
uncorrelated  with  each  other  and  uncorrelated  in  time.  Multiplying  this  equation  by  its  transpose  and 
taking  the  expectation  value,  one  finds 


Q(/  +  £|0  --  AAr. 


(39) 


Since  the  elements  of  the  state  vector  are  successive  integrals,  A  is  an  upper  right,  half  triangular 
matrix.  Equating  the  corresponding  elements  of  the  previous  equation,  one  finds 


A  (l  +  <5jf) 


(40) 


Sp'(t)fh 


0  yJSt{t)6  +  l\/SAt)S  +  f  y/s(  ( t)S 

0  0  y/s^)s  -  fv^ 

0  0  0  y/. St(t)S 
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The  algorithm  is  easily  mechanized.  A  random  number  generator  function  normally  produces  deviates 
uniformly  distributed  over  the  unit  interval.  They  may  be  converted  to  Normally  distributed  deviates 
using  one  of  the  standard  procedures.  Four  Normal  deviates  are  then  combined  to  form  the  noise 
vector  r  and  multiplied  by  A  to  produce  the  s  vector  needed  to  evaluate  Eq.  5  for  the  next  clock 
state.  The  element  u(t)  is  the  desired  simulated  clock  time  state.  The  simulated  clock  measurement 
is  obtained  by  adding  white  measurement  noise,  v(t),  to  u(t). 

z(t)  —  u(t)  +  v(t)  (41) 


TEST  RESULTS 

The  two  algorithms  were  run  simultaneously  in  order  to  intercompare  them.  Only  the  1987  version  of 
the  KAS-1  algorithm  (i.e.  without  automatic  parameter  estimation  and  mixed-mode  measurement 
capability)  was  available  for  testing. 

Figure  1  displays  the  difference  between  the  two  algorithms,  in  the  sense  ATI  -  KAS-1,  for  an  ensemble 
of  eleven  equally  weighted  clocks  with  no  data  rejection.  The  two  algorithms  are  seen  to  be  essentially 
identical  in  steady  state  on  the  long  term.  However,  as  mentioned  above,  ATI  is  not  designed  Tor  the 
natural  handling  of  transients,  while  KAS-1  is  still  optimal  in  such  situations. 

In  Figure  2,  ATI  was  allowed  to  exercise  its  3-standard-deviation  rejection  criterion,  while  the  KAS-1 
algorithm  was  forced  to  accept  all  the  data.  As  a  result,  on  the  short  term  ATI  is  shown  to  suffer 
from  small  discontinuities.  This  effect  is  avoided  with  the  KAS-1  algorithm  (in  its  usual  form)  by  its 
use  of  progressive  deweighting  of  outlying  data  based  on  Hampel’s  ip  functional. 

Figure  3  is  a  plot  of  Allan  deviation  vs.  sampling  time  t  for  one  of  the  eleven  clocks  (no  data  rejec¬ 
tion),  assuming  it  to  have  the  white  FM  and  random-walk  FM  noises  characteristic  of  a  conventional 
(HP5061A  model)  high-performance  cesium. 

Figure  4  is  a  similar  plot  for  the  ensemble  when  the  KAS— 1  algorithm  and  its  rejection  criterion  are 
employed.  Each  of  the  clocks  has  equal  white  FM  and  equal  random-walk  FM  noises;  hence,  they 
are  equally  weighted.  There  is  a  \/N  improvement  in  the  ensemble’s  collective  stability,  and  the  r  of 
minimum  Allan  variance  is  constant. 

Figures  5  and  6  show  simulated  data  of  contrasting  clock  types.  Figure  5  shows  the  frequency  stability 
of  VLGll-type  hydrogen  maser  without  a  long-term  frequency  drift  and  fig.  6  is  a  similar  plot  of  a 
SID  or  very  stable  (e.g.  a  De  Marchi-tuned)  cesium  (no  random  walk  of  frequency).  The  clock  types 
represented  in  figs.  3,  5,  and  6  are  used  in  an  ensemble  of  unequal  clocks  similar  to  that  which,  will 
soon  comprise  the  USNO  ensemble  (60  percent  cesiums,  30  percent  masers,  10  percent  STD’s),  This 
unequal  clock  set  is  used  as  the  input  to  two  algorithms  with  contrasting  weighting  philosophies. 

Figure  7  is  the  Allan  deviation  of  an  algorithm  (either  ATI  or  KAS-1)  which  automatically  determines 
the  weights  of  individual  clocks  in  the  ensemble  according  to  the  inverse  time  prediction  errors  at  the 
sampling  time;  whereas  fig.  8  is  the  deviation  of  an  algorithm  that  forces  all  clocks  to  be  weighted 
equally.  Comparing  the  two  figures,  representing  eight  years  of  data,  the  weighted  algorithm  is  signif¬ 
icantly  better  in  the  short-term,  whereas  the  equally  weighted  algorithm  is  better  in  the  long-term. 
Future  versions  of  the  KAS  algorithm  will  produce  optimum  results  at  all  sampling  times  by  completely 
accounting  for  the  clock-ensemble  correlations. 

Figures  9  and  10  are  analogous  to  Figs.  7  and  8,  respectively,  but  display  25  years  of  data.  The  same 
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effects  are  noted,  with  the  long-term  degradation  of  the  weighted  algorithm  being  even  more  evident. 
Clearly,  neither  the  ATI  nor  the  KAS-1  algorithm  when  using  a  weighting  scheme  is  optimal  for  sets 
of  dissimilar  clocks.  A  new  algorithm,  KAS-2,  has  been  developed  by  Ball  Aerospace  that  is  optimum 
for  all  clock  types  and  all  sample  times.  It  will  be  tested  in  the  Algorithm.  Test  Bed  in  1990. 

CONCLUSION 

The  use  of  simulated  data  is  shown  to  be  a  powerful  method  of  comparing  different  clock  models  and 
ensembling  algorithms.  A  Kalman  filter-based  algorithm,  such  as  KAS  -1,  can  use  unequally  spaced 
data  and  can  handle  transients  well,  unlike  an  ad  hoc  procedure  like  the  ATI  (1987)  algorithm.  The 
KAS-1  algorithm  is  able  to  handle  a  mix  of  time  and  frequency  measurements  quite  naturally. 

Data  rejection  with  ATI  leads  to  small  discontinuities  in  the  timescale,  unlike  the  case  for  KAS-1. 
While  both  algorithms  are  optimal  in  steady  state  for  clocks  having  proportional  noises,  they  are  not 
optimal  for  sets  of  significantly  different  clocks.  The  stability  of  a  real  clock  ensemble  would  be  further 
degraded  by  clock  intercorrelations  and  systematic  errors,  so  that,  if  the  clock  weights  are  based  only 
on  Allan  variances,  the  results  may  not  be  as  good  as  those  employed  simply  using  equal  weights^]. 

As  the  project  continues,  more  study  will  be  made  of  simulated  data;  real  clock  data  will  be  analyzed; 
and  the  improvements  to  the  KAS-1  algorithm  suggested  herein  will  be  implemented  and  tested. 
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Figure  1:  ATI  and  KAS-1  are  essentially  identical  in  steady  statf 
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Figure  2:  Small  discontinuities  in  ATI  rejection  method. 
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Figure  3:  Simulated  cesium  clock  data. 
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Figure  4:  Square-root  of  N  improvement  with  11  equal  clocks. 
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Figure  7:  Output  of  algorithm  using  inverse  variance  weighting  (8  years  of 


simulated  data). 
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Figure  8:  Output  of  algorithm  using  equal  weights  for  all  clocks  (8  years  of  simulated  data) 
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Figure  9:  Output  of  algorithm  using  inverse  variance  weighting  (25  years  of  sin 
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Figure  10:  Output  of  algorithm  using  equal  weights  for  all  clocks  (25  years  of  simulated  data) 
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